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AN IRREVERSIBLE CONSTITUTIVE LAW FOR MODELING THE DELAMINATION 
PROCESS USING INTERFACE ELEMENTS 

Vinay K. Goyal* * * § and Eric R. Johnson* 

Virginia Polytechnic Institute and State University, Blacksburg, VA 24061-0203 
Carlos G. Davila* 

NASA Langley Research Center, Hampton, VA, 23681 

Navin Jaunky 5 

Institute for Computer Applications in Science and Engineering, Hampton, VA, 23681 

An irreversible constitutive law is postulated for the formulation of interface elements to predict initiation 
and progression of delamination in composite structures. An exponential function is used for the constitutive 
law such that it satisfies a multi-axial stress criterion for the onset of delamination, and satisfies a mixed mode 
fracture criterion for the progression of delamination. A damage parameter is included to prevent the restoration 
of the previous cohesive state between the interfacial surfaces. To demonstrate the irreversibility capability of 
the constitutive law, steady-state crack growth is simulated for quasi-static loading-unloading cycle of various 
fracture test specimens. 


INTRODUCTION 

Delamination in composite structures usually originates 
from geometric discontinuities and material defects such as 
free edges, dropped plies, re-entrant comers, notches, and 
transverse matrix cracks. Recently, significant progress has 
been made in the development of tools to predict intralam- 
inar damage, which often precedes the onset of delamina- 
tion. Delamination can be a major failure mode in compos- 
ites structures and can lead to significant loss of structural 
integrity. The virtual crack closure technique (VCCT ) 1,2 
has been successfully used in the prediction of delamination 
growth. However, an initial delaminated area must be prede- 
fined and a self-similar crack growth is assumed. 

To overcome the limitations associated with the VCCT, 
interface elements can be located between composite lam- 
ina to simulate initiation of delamination and non-self-similar 
growth of delamination cracks without specifying an initial 
crack. Delamination is initiated when the interlaminar trac- 
tion attains the maximum interfacial strength, and the de- 
lamination front is advanced when the local surface fracture 
energy is consumed. A softening constitutive law that relates 
tractions to the relative displacements is generally used to for- 
mulate interface elements. The softening constitutive law is 
based on the Dugdale 3 and Barenblatt 4 cohesive zone model 
to expunge the singular stress field ahead of the crack-tip en- 
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countered in linear elastic fracture mechanics. The softening 
portion of the constitutive law models the degradation of the 
material ahead of the crack-tip. For laminated composites 
this degradation includes nucleation, growth and coalescence 
of microcavities. Hilleborg 5 developed the first comprehen- 
sive interface finite element model and applied this method in 
concrete cracking. Later, Needleman 6 developed a cohesive- 
decohesive formulation to simulate dynamic crack growth in 
isotropic elastic solids. 

The exact mathematical form of the interfacial constitu- 
tive law is less important than its capability to represent the 
maximum interfacial strength and critical fracture energy. 
Functions with continuous derivatives have a numerical ad- 
vantage over functions with discontinuous derivatives when 
used with Newton-Raphson method because the tangent stiff- 
ness is smooth. A smooth tangent stiffness as a function of 
the relative opening displacement has been found to mitigate 
the numerical oscillations encountered in using a softening 
constitutive relation and to eliminate oscillatory convergence 
difficulties 7 . 

The exponential function for the softening constitutive law 
is smooth and mimics the physics involved in the separation 
of two atoms initially bonded 3 . This form of the constitutive 
law has been used in the analysis of crack initiation, dynamic 
growth, branching, and arrest in homogeneous materials 9 . 
Shahwan and Waas 10 used it to study delamination of com- 
posite structures caused by compressive loads. The various 
exponential constitutive laws that have been successfully em- 
ployed to simulate delamination are based on the assumption 
that the consumed local surface fracture energy can be recov- 
ered. This assumption is not valid for structural systems with 
stresses that may internally redistribute upon external load- 
ing. The cracks may arrest and cracks faces may close. Ortiz 
and Pandolfi 11 postulated a damage model and used an expo- 
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nential constitutive law to account for such irreversibilities. 
A limitation of this model is that the critical energy release 
rates and the maximum interfacial strengths associated with 
Mode I, Mode II, and Mode III fracture cannot be specified 
separately. 

The present work aims at the establishment of an exponen- 
tial softening constitutive law that satisfies empirical mixed- 
mode delamination failure criteria for the onset and progres- 
sion of delamination. An internal state variable is included 
in the constitutive law to permanently damage the internal 
surfaces that have exceeded maximum strength during the de- 
formation process. The paper is structured as follows: (i) 
mixed-mode fracture criteria, (i) mechanics of interfacial sur- 
faces, (ii) interface finite element, (iii) finite element results, 
and (iv) concluding remarks. 

MIXED-MODE FAILURE CRITERIA 

A quadratic failure criterion based on interlaminar tractions 
has been used to predict onset of delamination 12 . To simu- 
late the progression of delamination under mixed-mode load- 
ing conditions, the power law form of the fracture criterion 
that includes Mode I, Mode II and Mode III interaction has 
been successfully used with a bilinear constitutive law 13-15 . 
Davila and Camanho 16 developed a bilinear constitutive law 
that can be used with any mixed-mode failure criterion 16 . To 
the authors’ knowledge, no work has been found incorporat- 
ing empirical failure criteria into the exponential softening 
constitutive law. A brief description of the failure criteria 
used in this paper is presented next. 


Criterion for Progression of Delamination 

Delamination propagates when the energy release rate 
equals its critical value under pure Mode I, Mode II, or 
pure Mode III fracture. Generally, delamination growth oc- 
curs under mixed-mode loading. Under this type of load- 
ing, delamination growth might occur before any of the en- 
ergy release rate components attains its individual critical 
value. The power law criterion based on the one proposed 
by Whitcomb 18 is 




+ 



1 (2) 


where Gj is the energy release rate under Mode j fracture, 
and Gj c is the single-mode critical energy release rate for 
j = I, II, III. The material parameter a defines the shape 
of the failure locus. For a = 1, one recovers the linear 
interaction criterion 19 . The shape of the failure locus is a tri- 
angular surface. The shape of the failure surface approaches 
a 1 /8-cube surface as a increases from 2. Reeder 20 evaluated 
different fracture criteria for mixed-mode delamination in a 
brittle graphite/epoxy composite, a toughened graphite/epoxy 
composite, and a tough graphite/thermoplastic composite us- 
ing the mixed-mode bending (MMB) test specimen. The 
power law criterion was a reasonable fit to the test data for 
the three different materials. Thus, the failure criterion in 
Equation (2) is incorporated into the constitutive law of the 
interface material. 


MECHANICS OF THE INTEKKACIAE SURFACES 


Criterion for the Onset of Delamination 

Under pure Mode I, Mode II, or pure Mode III loading, the 
onset of delamination occurs when the corresponding inter- 
laminar traction exceeds its respective maximum interfacial 
strength. However, under mixed-mode loading, delamination 
onset may occur before each traction component reaches its 
maximum interfacial strength. An expression that considers 
the interaction between the traction components under mixed- 
mode loading is the multi-axial stress criterion given as 



where Tj is the interlaminar traction component associated 
with the j-direction, Tj is the maximum interlaminar trac- 
tion, and {£) = £, if t > 0, otherwise it is zero. This function 
has been included to emphasize that the normal compressive 
traction X3 does not contribute to the onset of delamination. 
In Equation 2 ,T e is an effective normalized traction, and 
a > 2 is a real number that determines the shape of the 
tri-dimensional failure surface. The quadratic delamination 
interaction is recovered from Equation (1) with a = 2. The 
failure surface for a = 2 is a convex semi-sphere in the space 
of normalized tractions T.JTf j = 1, 2, 3. As the value of 
a is increased, the failure surface approaches a half-cube sur- 
face. 


Interfacial surfaces consists of an upper surface S + and 
lower surface S . The upper surface corresponds to the up- 
per bulk material, and the lower surface corresponds to the 
lower bulk material. The surfaces S ± are coincident with a 
reference surface S° in the undeformed configuration as is 
shown in Figure 1 . Thus, it is said that the interface material 
is of zero thickness. The surfaces S ± independently displace 
and stretch, and are connected by a continuous distribution 
of nonlinear springs that act to resist the Mode I opening or 
Mode II and Mode III sliding of the upper and lower surface. 

It is convenient to define a mid-surface S m where the trac- 
tions and relative displacements are evaluated. For this pur- 
pose, let us consider any two points P + and P~ contained 
in S + and S~ and coincident in the undeformed configura- 
tion. The locus of the midpoints P m of the line joining P + 
and P ~ define the mid-surface S m of the interface material. 
Refer to Figure 2. The normal and tangential components of 
the traction and relative displacement vector are determined 
by the local orientation of the mid-surface S m . The virtual 
work done by the cohesive-decohesive tractions is given by 

5Wiut = JJ s SAj Tj dS m (3) 

for any kinematically admissible relative displacements A j, 
where Tj are the interlaminar traction components acting on 
a unit deformed area conjugate to the relative displacements, 


2 OF 12 


American Institute of Aeronautics and Astronautics Paper 2002-1576 




Fig. 1 Interface material deformation. 



*1 


Fig. 2 Interface material mid-surface. 

and S m is the surface area. The resistive tractions that are 
associated to the relative displacements at the point P m are 
shown in Figure 2. The interlaminar normal traction is de- 
noted T 3 and the tangential tractions are denoted 7\ and T 2 . 

In the next section, the components of the relative displace- 
ments are obtained in terms of the displacement field with 
respect to the undeformed configuration. Next, the consti- 
tutive equations that relate the relative displacements to the 
traction field are presented. The kinematics and the con- 
stitutive modeling fully describe the mechanics of interface 
debonding. 

Kinematics of the Interface Material 

The fundamental problem introduced by the interface ma- 
terial is the question of how to express the virtual relative 
displacements between the surfaces S ± in terms of virtual 
displacements. As shown in Figure 1, consider a three- 
dimensional space with Cartesian coordinates X t , * = 1,2, 3, 
and let there be surfaces S ± coincident with S° defined in 
this space by X ; = X. t (i]i, 7/2 ) , where 771,772 are curvilinear 
coordinates on the surface S' 0 . 

Let the Cartesian coordinates xf = xf (771 , 772), * = 1, 2, 3 
describe motion of the upper and lower surfaces S ± in the 
deformed configuration. Any point on S ± in the deformed 
configuration is related to the same point on S° through 

xf = Xi + Uf (4) 


where Uf are displacement quantities with respect to the 
fixed Cartesian coordinate system. The coordinates x™ = 
x? (771,772),* = 1,2,3 define the mid-surface S m given by 

x? = \ (xf + xt) = X, + \(Uf + Ur) (5) 

The surface S m is coincident with S° in the undeformed 
configuration. As mentioned earlier, the components of the 
relative displacement vector are evaluated at the mid-surface 
S m . Therefore, the local orientation of normal and tangential 
unit vectors to the surface S m is required. This is, 


ri 


dxf dxf 1 T 
din ’ 9771 ’ din J 


( 6 ) 


, _ (chef dxf dxf\ T 
1-2 \ dr ] 2 ’ 0?72 ’ 0772 J 

and the normal vector is simply 


(7) 


r 3 = ri x r' 2 (8) 

The tangential vectors r 3 , r' 2 may not be perpendicular in a 
curvilinear coordinate system so that, 


r 2 = r 3 x ri 


(9) 


For i = 1, 2, 3, the normal and tangential unit vectors to the 
surface S m at a point P m € S m are 


hi = 



( 10 ) 


These unit vectors define the local orthogonal coordinate 
system at S m and is related to the fixed coordinate system 
through the rotation matrix 


R=[ri,r 2 ,r 3 ] (11) 

The normal and tangential components of the relative dis- 
placement vector expressed in terms of the displacement field 
is, 

A» = Rjiixj - xj) = Rji (Uf - Ut) (12) 

where Rji are components of the rotation matrix. Since x™ 
depends on the displacements Uf, the rotation matrix also 
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depends on Uf. Therefore, the virtual relative displacement 
are expressed in terms of the virtual displacements as follows. 


5A,; =(% + [/+ 


SUf- 


. dRki 


Rh + U, 

3 k dU- 


SU 7 


5A i = Q±5U+-Qj l 5Ur (13) 

Equation (13) is substituted into Equation (3) to obtain the 
expression of the internal virtual work in terms of the virtual 
displacements. This form of the internal virtual work is con- 
venient for the finite element formulation. In addition, the 
differential surface area of the mid-surface dS m in the de- 
formed configuration is expressed in the form, 

dS m = MdS° (14) 


where M is a function of the displacement field Uf, and dS° 
is the differential undeformed surface area. 


Constitutive Equations for the Interface Material 



Fig. 3 Traction-Stretching curve of spring as a function of p 

a high spring stiffness maintains the points together. Under 
isothermal conditions, the traction T that acts to resist the 
stretching A of the spring is expressed as 


The stress singularities at the crack-tip in the linear elas- 
ticity solutions, stemming from the sharp slit approximation, 
cannot be reconciled with any realistic local rupture pro- 
cess. From the molecular theory of strength it is known that 
there exists stress limits for which molecular bond rupture 
occurs. The softening-type of cohesive zone model is in- 
tended to represent the degradation of the material ahead of 
the crack-tip. It captures strength-based bond weakening, and 
fracture-based bond rupture. The mechanics of the delam- 
ination process comprises three interrelated phases: (i) the 
initiation of delamination, (ii) the evolution of the degrada- 
tion zone, (iii) and the delamination growth. The first phase 
that takes place is the initiation of delamination, and it is 
based on a stress limit determined experimentally. A stress 
measure that is used as the limiting value, may involve an in- 
teraction of interlaminar stresses such as the equivalent Von 
Mises stress, or that in Equation (1). The second event is the 
development of a zone ahead of the crack-tip that experiences 
intense deformation such as plastic deformation in metals, 
elongated voids that contains a fibrous structure bridging the 
crack faces in polymers (crazing), and high density of tiny 
cracks in brittle ceramics. The molecular bonds are weakened 
and the nonlinear softening behavior is confined in this degra- 
dation zone, or process zone. The third event, is the growth 
of delamination, bond-rupture, and it is based on a fracture 
criteria such as Equation (2). The constitutive equations to be 
developed, mathematically describe these three delamination 
phases. The focus of this section is to develop the consti- 
tutive equation for single-bond rupture based on continuum 
damage mechanics approach. This particular case is extended 
to mixed-mode delamination. The constitutive equations that 
are postulated in this section, are shown to satisfy the failure 
criteria for initiation and progression of delamination pre- 
sented in the previous section. 

Let assume that the two points P + and P~ contained in 
S + and S~ as shown in Figure 2 are connected with a spring. 
The points are coincident when the spring is unstretched, and 


T( A) = T c A exp 


1 - A 1 3 

~T~ 


(15) 


where A = A/A°, and T c is the maximum bonding strength 
that occurs at the critical stretching value A 0 . The parameter 
6 with p > 1 and p £ R + defines the stretching range for 
which the bond is weakened before complete rupture occurs. 
It is in this range, that damage accumulates. In Figure 3, the 
traction-stretching curve is shown for different values of the 
parameter 0. The work of debonding per unit area, G c , is 
given by the area under the traction-stretching curve, or, 


G c 


f 


T(A)dA 


T c A c pv-pnv r 


2 

P 



(16) 


Tfy] is the Euler gamma function of 2 , and F[l/2] = ff. 
By prescribing T c , G c , and p in Equation (16), the parameter 
A c can be computed. The exponential function in Equation 
( 1 5) is a suitable representation of a softening constitutive law 
because with increasing stretching of the spring A, the trac- 
tion T increases to a peak value T c and then decreases until 
complete debonding occurs. Equation (15) is only valid for 
monotonically increasing separation because the consumed 
debonding energy can be recovered upon unloading. 

An internal state variable d that tracks the damage state of 
the spring needs to be included in Equation (15) to account 
for irreversible effects. In the following irreversible law an 
elastic damage model instead of a plastic damage model is 
assumed, 


T( A) =T C A exp 


2- A 13 /d-d 

P 


(17) 


Within the framework of continuum damage mechanics, it 
is possible to impose restrictions on d. It must increase as a 
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Fig. 4 Traction-stretching curve as a function of the evolution 
of damage of the spring with (3 = 1 


In reference to Figure 2, the components of the normalized 
relative displacements between P ± with respect to the orien- 
tation of the surface S m at a point P m is, 

v = Ai ii + A 2 i 2 + A 3 i 3 (20) 

where ii, i 2 , 13 are the unit vectors normal and tangent to the 
surface S m at a point P rn . An effective relative displacement 
A is defined by the norm of v 

A = y / A^~+~A|~+ _ A^ ( 21 ) 

We assume that the normalized scalar traction T v acts along 
the direction of v to resist the effective relative displacement 
A . The proposed constitutive law for the interface material is 
defined along v, 

T v (Ai, A 2 , A 3 ) = A Q(Ai, A 2 , A 3 ) ( 22 ) 


function of time because thermodynamics requires that the ir- 
reversible dissipation associated with the debonding process 
remains semi-positive, i.e., d > 0. An equivalent mathemati- 
cal expression is 

max^M^-^A^) , d (0) = 1 (18) 

with ti > ti_i. If the spring is assumed undamaged at t = 
to, then the initial condition is d ^ = 1. Equation (17) is 
equivalent to Equation (15) if no damage occurs, d = 1, or 
for monotonic increasing loading, d = AA Unloading does 
not occur linearly to the origin, but with an exponential form. 
The energy of dissipation associated to fatigue is neglected 
in this work. This assumption is valid in the case of a spring 
that undergoes a small number of loading-unloading cycles. 
Thus, future work will be aimed at extending the Equation 
( 18) to incorporate fatigue. 

Equations (17) and (18) with (3 = 1 are used for the 
traction-stretching curve in Figure 4. The labels 1, ..., 6 in this 
figure, represent the damage evolution of the spring connect- 
ing P ± . The spring is unstretched at point 1 . With increasing 
stretching, a cohesive traction develops to resist the separa- 
tion. At point 2, the spring stiffness holds P ± together in the 
quasi-linear range of the law. The onset of delamination oc- 
curs at point 3, where the traction attains its maximum value. 
As the spring is stretched beyond the onset of delamination 
to point 4, damage is accumulated in the spring and the trac- 
tion gradually decreases. The spring is partially unstretched 
from point 4 to point 5, and unloading occurs. The spring is 
stretched again to point 6, and the loading traction-separation 
curve is exactly retraced upon unloading. The traction even- 
tually vanishes as the spring is stretched. 

Equations (17) and (18) are extended to the mixed-mode 
delamination case. To develop the constitutive equations, it 
is convenient to normalize the relative displacements A j and 
the tractions Tj with respect to the critical separation values 
Aj and the maximum interfacial strengths T?, 

A; A,-/A7, 7) (19) 


where Q is a decreasing function of any of the normalized 
relative displacements A j, j = 1,2,3. The components of 
the traction acting along v, normal and tangent to the mid- 
surface S m at a point P ,n is 

Tj = T v i j ■ jj^| = A j Q (23) 

for j = 1,2,3. The function Q is chosen to satisfy the 
multi-axial stress criterion in Equation 1 for the onset of de- 
lamination and the mixed-mode fracture criterion in Equation 
2 and is given by 


Q = exp 


2-^ /d-d 

0 


(24) 


with a scalar mixed-mode parameter /z that couples the nor- 
malized relative displacements for the opening and sliding 
mode 

P = (|A 1 |“ + |A 2 | a + {A 3 )“) 1/ “ (25) 

where |-| is the absolute value function, and {£) = £ if £ > 0, 
otherwise it is zero. The material parameter a defines the 
shape of the failure surface for the onset and progression of 
delamination. The internal state variable d is given by, 

//(u) _ max ^ fjPq'j , d*- 0 ) = 1 (26) 


The constitutive equations are slightly modified to take into 
consideration the mechanical behavior of the interface mate- 
rial under contact conditions. The surfaces S± are assumed 
smooth so that frictional effects can be neglected. When con- 
tact is formed between two smooth surfaces, the equilibrium 
largely depends upon the distribution of elastic forces in the 
contacting surfaces. Two surfaces are under contact at a point 
pm , pm g s m if the relative displacement A 3 between P ± 
is less than zero. For A 3 < 0, a large repulsive traction T 3 
develops to avoid interpenetration of the surfaces S ± at P m . 
The constitutive equations for mixed-mode delamination are 
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obtained from Equations (23) to (26), and summarized as fol- 
lows 


A 2 = £3 A 3 with | 2 and £3 fixed during the loading history. 
The terms in Equation (2) are evaluated as follows, 


( 


Ti 

T 2 

t 3 


+ 


Ai 

A2 

(As) 



2 - p? /d-d 

P 


(27) 


0 

0 


-(-As) 




( Jo 3 A?(Ai, A 2 , A 3 )dA 3 \ 
V fo°°T 3 (0,0,A 3 )dA 3 J 


1 

(!+£? + $ 


\2/c 


+ </>l(A 3 


*/2 


and k,k> 1 is an interpenetration factor to magnify the re- 
pulsive force T 3 , and chosen arbitrarily. Equations (26) and 
(27) reduce to Equations (17) and (18) for single-mode de- 
lamination. 

The empirical parameters governing the constitutive 
equations in (27) are the critical energy release rates G/ c , 
Giic, Gnic', the maximum interfacial strengths Tf, T 2 , T3 ; 
and the the critical separation values Af , A 2 , Ajj. These may 
be specified based on atomistic models of separation or on 
a phenomenological basis depending whether the separation 
process is governed by ductile void coalescence or a brittle 
cleavage mechanism. By specifying the critical energy 
release rates and the maximum interfacial strengths, one can 
obtain the critical separation values. The path independent J- 
integral along a boundary that contains the interface material 
can be used to show that the area under the traction versus 
separation curve is the work of fracture per unit area. Equa- 
tion (16) under pure Mode I, Mode II, or Mode III fracture, 
is used to obtain the critical separation values A'f j = 1, 2, 3. 



// 0 Al T 1 (A 1 ,A 2 ,A 3 )dA 1 \ a/2 

V / 0 oo Ti(A 1 , 0 , 0 )dA 1 J 


& 

(i + £? + £§ 


,2/a 


+ (j>2 (A 3 


*/2 



/ f 0 A2 T 2 (A 1 ,A 2 ,A 3 )dA 3 \ a/2 
V fo° T 2 (0, A 2 , 0 )dA 2 J 




(1 + ?2 a + tf) 


y^ + MAa) 


a/2 


where <j>j(A 3 ),j = 1 , 2, 3 are exponential decaying functions 
with increasing A 3 . The progression of delamination occurs 
when the functions <j>fy A 3 ) , j = 1 , 2 , 3 are virtually zero. 
Adding the last three equations shows that the power criterion 
in Equation (2) is satisfied. ■ 


Proof. The exponential constitutive law in Equations (26) 
and (27) satisfy Equation (1) for the onset of delamination, 
and Equation (2) for the progression of delamination. 

For simplicity, monotonically increasing loading is as- 
sumed, i.e., d = pP . The effect of interpenetration is also 
neglected, A 3 > 0. For the onset of delamination, the com- 
ponents of the traction vector in Equation (27) are substituted 
into Equation (1) to obtain the effective traction T e , 



This equation is analogous to Equation (15) for single-mode 
delamination. In view of Equation (28), delamination onset 
occurs when p = 1. At this value of p, the effective traction 
attains the maximum value of one. The failure criterion in 
Equation (1) predicts delamination onset at an effective trac- 
tion equal to one. Therefore, with the proposed constitutive 
law in Equation (27) delamination initiates when the criterion 
in Equation (1) is satisfied. 

For the progression of delamination, proportional straining 
is assumed. The relative displacement associated to the slid- 
ing Mode II and Mode III are written as Ai = £ 2 A 3 and 


IN TERFACE FINITE ELEMENT 

The formulation for the interface element is based on the 
work of Beer 21 . A non-linear solution procedure is necessary 
because of the geometrical nonlinearities and the nonlinear 
mechanical behavior of the interface material. The objective 
of this section is to obtain the tangent stiffness matrix K® and 
the internal force vector f? t required in the nonlinear solution 
procedure. 

A 2n-noded isoparametric interface element with 6 ?? de- 
grees of freedom and applicable to three-dimensional analysis 
is used. The element consists of an upper and lower surface 
Sf with n-nodes each. The natural coordinate system is rp 
and rj 2 . For the surfaces Sp, node j has three translational 
degrees of freedom qf-, q 2 j,q 3 j with the first subscript imply- 
ing the associated global direction. The nodal displacement 
vector q is arranged as follows, 

q={q + ,q~} T (29) 

q = {■■aqipqfpqfj ,-} 7 

and j denotes the node number, j = 1, ...,n. The displace- 
ment field UP(t) i,r/ 2 ), j = 1,2,3 for the surfaces Sf are 
independent and in terms of the global displacement degrees 
of freedom qf- 

uf(m,V2) = qfn N n(q urn) (30) 
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where N n is the shape function corresponding to the n - th de- 
gree of freedom. Substituting Equation (30) into Equation 
(13) gives. 


Material Tangent Stiffness 

The components of the material tangent stiffness D are oh 
tained in the incremental form, 


5Ai = QjiN n 5qf n - QjiN n 5q in ( 31 ) 

Equation (31) in matrix form is, 

5 A = [C£ N, -Q t N] 5q (32) 

= [B+,-B ] 5q= Bdq 

where N is 

N = [..., Nj I, ...], j = 1, ...,n (33) 

and I is a 3 x 3 identity matrix. Equation (32) relates the 
relative displacement to the nodal displacement degrees of 
freedom. 

The internal force vector of the interface element is ob- 
tained by substituting Equation (32) in (3), 


5W? nt = Sq T [[ B t T dS. r = d'q T f^ t (34) 
J Js™ 

where T is the traction vector acting on the deformed mid- 
surface and the integration is performed over the deformed 
element mid-surface. In numerical analyses, the internal 
force vector needs to be computed accurately, and the tan- 
gent stiffness matrix may be computed approximately. The 
computation of the tangent stiffness matrix is intensive and 
a very accurate expression is not required. Therefore, the 
partial derivatives of the differential area in Equation (14) is 
neglected. For the computation of KJr, the derivatives of the 
rotation matrix with respect to the nodal displacements are 
neglected. This approximation with Equation (32) leads to 


B+ = B- = B s (35) 

5 A = [B s , -B s ] d'q = B'<5q 

Thus, the approximate tangent stiffness matrix is, 

Kf = — ^ « [[ B ,T DB ' dS™ (36) 

C>q JJs T 


where D is the material tangent stiffness, and is later defined. 
Equation (36) is rewritten using the relation in Equation (35), 


where 


K 


K s — K s 
— K s K s 


K s 



B s t DB s dS™ 


(37) 


(38) 


The internal force vector is accurately computed, while the 
approximations for the tangent stiffness matrix save compu- 
tational time because only a quarter of the full matrix has to 
be computed. 


(7 / • 

5Ti = g^SAj = Dj.jSAj (39) 

First consider the case in which there is no interpenetration, 
that is, for A3 > 0. The components of D are obtained by 
differentiation of Equation (27) according to Equation (39), 


Da 



^%iA 3 r 2 


Q 


(40) 


where % is the Kronecker delta, Q is given by Equation (24), 
and ui is defined by, 


Co 


1 if d = iff 
d if d> ^ 


(41) 


Now consider the case for which interpenetration is de- 
tected, that is, A3 < 0. The non-zero components of D are 
given by Equation (40) for i,j = 1,2 and the component re- 
lated to interpenetration, 

D 33 = K 0 (l + k [Aal' 3 ) exp ^ ^ (42) 

where K 0 = T§ exp(l//?)/Ajj. The range of the values 
of D 33 should be restricted by two conditions: (1) A small 
D 33 induces interpenetration, and (2) a large D 33 produces 
ill-conditioned matrices. A list of references on these restric- 
tions is given by Davila et al. 22 . The value of D 33 should be 
in the range. 


10 6 N/mm 3 < D 33 < 10 7 N/mm 3 

The upper bound of the condition cannot be guaranteed be- 
cause of the exponential nature of Equation (42). Therefore, 
for A3 < 0, the expressions T 3 and D 33 are modified to have 
the form 

T 3 = K 0 A 3 , D 33 = Kq (43) 

and A'o = T| exp(l/(3)/A$. 

The material tangent stiffness is non-symmetric, and can be 
positive definite, semi-definite, or negative definite. For ji > 
1, the matrix D t j is negative definite. The material tangent 
stiffness matrix has properties of an anisotropic material, one 
which has strong dependence on the relative displacements 
in all directions. For single-mode delamination, D is fully 
diagonal, otherwise, some of the off-diagonals are non-zero. 

Consistent and Inconsistent Tangent Stiffness 

For the full-Newton-Raphson nonlinear solution proce- 
dure, the consistent tangent stiffness matrix is used in the 
finite element analysis. However, when softening constitu- 
tive laws with the consistent tangent stiffness are employed, 
the tangent stiffness matrix is often ill-conditioned and a 
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converged solution may not be obtained 23 . An alternative 
solution is to refine the mesh ahead of the crack-tip or the 
decrease the maximum interfacial strength 14, 24 . Refining the 
mesh size increases the computational time, and lowering the 
maximum interface strength can result in a premature initia- 
tion of delamination 7 . Alternatively, researchers often utilize 
a positive definite matrix such as the material secant stiffness 
when dealing with softening constitutive laws. However, a 
large number of iterations results in using the material secant 
stiffness. As an alternative, three different modifications to 
the tangent stiffness matrix eliminate these convergence diffi- 
culties while a converged solution can be obtained in a small 
number of iterations: 

1. Equation (36), Kf ; = max(0, Kf { ), i = 1,2, ..., 2 n 

2. Equation (37), K Sii = max(0, K Sii ), i = 1,2, ..., n 

3. Equation (40), Du = max(0, Du), i = 1, 2, 3 

The convergence rate of option 1 is better than option 2, and 
the convergence rate of option 2 is better than option 3. If the 
mesh is coarse, is better to choose option 3. 

Contact Elements 

Interface elements were developed to model initial delam- 
inated surfaces. All the components of the material tangent 
stiffness is zero, except for the case in which interpenetra- 
tion is detected. If interpenetration is detected Equation (43) 
is used. Thus, these interface elements act like contact ele- 
ments. 

FINITE ELEMENT RESULT S 

Numerical results are presented for quasi-static loading and 
unloading of the double cantilever beam (DCB), the end load 
split (ELS), end notch flexure (ENF), and fixed ratio mixed 
mode (FRMM) fracture test specimens. Results are also pre- 
sented for quasi-static loading of the mixed mode bending 
(MMB). Mode I fracture occurs in the DCB specimen. Mode 
II occurs in the ELS and ENF specimens, and Mode I and II 
occur in the FRMM and MMB. The fracture test specimens 
are shown in Figure 5. 

Mode I and mixed-mode test specimens are modeled with 
the laminate stacking sequence [Of!] and the unidirectional 
material properties of Graphite -Epoxy listed in Table 1 . An 
isotropic material with E = E n and v = v \2 are used for the 
Mode II test specimens rather than composite. The maximum 
interfacial strength and the critical energy release rates are 
listed in Table 2. The geometrical properties are the length 
L = 100 mm, the arm thickness h = 1.5 mm, and width 
B = 10 mm. For the DCB, the geometrical properties are dif- 
ferent from the other test specimens: L = 150 mm, h = 1.5 
mm, and B = 20 mm. The initial crack length ao of each test 
specimen is: DCB - 50 mm, ENF - 30 mm, ELS - 50 mm, 
FRMM - 40 mm, and MMB - 20 mm. 

The interface elements are positioned between the up- 
per 0° laminate and the lower 0° laminate. Delamination 
is constrained to grow in the plane between the upper and 


Table 1 Graphite-Epoxy Properties 


E n 

E 22 , E 

33 C 12 , G 13 

Q) 

< 

11 

< 

>5 

150.0 GPa 

11.0 GPa 6.0 GPa 

3.7 GPa 0.25 

0.45 


Table 2 

Interface Material Properties 


Ti, T 2 

D 

Gic 

^IIC 5 GfflC 

K h 

80 MPa 

60 MPa 

0.352 N/mm 

1 .45 N/mm 10 7 

N/mm 3 


lower laminates. Interface elements with contact properties 
were placed along the initial crack length and interface el- 
ements formulated with the softening law are placed along 
the bonded length. The upper and lower laminates are mod- 
eled with C3D8I incompatible-mode 8 node solid element 
available in ABAQUS. Each laminate is modeled with one 
element through the thickness, 100 elements along the length 
of the laminate, and one element across the width. See Fig- 
ure 6a. For the DCB, three elements along the width are 
used. The eight node isoparametric interface element for 
three-dimensional analysis shown in Figure 6b is compatible 
with C3D8I solid element. The element was implemented 
in the commercial finite element code ABAQUS as an UEL 
subroutine. Three point Gauss integration is used for the 
computation of the tangent stiffness matrix and internal force 
vector. 

An incremental-iterative approach is adopted for the non- 
linear finite element analysis, and the Newton’s method avail- 
able in ABAQUS is used to trace the loading path of the spec- 
imens with a displacement-control analysis. For the MMB, 
the Riks method available in ABAQUS is used. The modifi- 
cation to the tangent stiffness matrix mostly used is option 2 
discussed in the section of interface elements. The response 
of the test specimens is characterized by the load-deflection 
curve. A typical finite element model of one of the test spec- 
imens consists of about 300 elements, and 2000 degrees of 
freedom. The computational time required was about 1200 
seconds of CPU time on a Sun Solaris 2000. The average 
number of iterations per load increment is 7. 

The finite element solutions are compared to the beam 
analytical solutions derived from linear elastic fracture me- 
chanics. The analytical solutions for the DCB and ENF are 
given by Mi et al. 14 , and for the FRMM and ELS are given 
by Chen et al. 24 . The finite element solutions for the MMB 
test specimen are compared to the analytical solution in the 
appendix. 

The DCB test specimen shown in Figure 5a is used to de- 
termine the interlaminar fracture toughness in Mode I. The 
load w is symmetrically applied, equal and opposite at the 
tip of the upper and lower arm of the DCB test specimen. 
The corresponding reaction force P is computed. The other 
end of the specimen is clamped. The response of the DCB is 


8 OF 12 


American Institute of Aeronautics and Astronautics Paper 2002-1576 



P, w 



Pi, w 2 P\, Wi 

^ |< Cld S’ 




(c) Mode II - ENF, P, = 0 

Mixed-Mode - MMB, P 2 t- 0 
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(b) Mode II - ELS 


N d!o ■ 


P, W 


P, w 

A 


(d) Mixed-Mode - FRMM 

Fig. 5 Fracture test specimens. 
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laminate 



Interface Lower 

elements laminate 


a) Finite element model of the ENF test specimen 
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b) Eight-node isoparametric interface element 


Fig. 6 Finite element modeling 
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Direction of edge 


crack growth 


a) Load-opening response of the DCB 


b) Non-self similar delamination growth 


Fig. 7 DCB specimen with ao = 50mm 
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Fig. 8 Load-displacement response of Mode II test specimens. 
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Fig. 9 Load-displacement response of the FRMM test specimens 
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Fig. 10 Load-displacement response of the MMB test specimens 
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shown in Figure 7a. For a loading-unloading cycle, excellent 
agreement of the FEM results are obtained compared to to 
the closed form solutions and to the experimental data. A top 
view of the Mode I specimen near the delamination front is 
shown in Figure 7b. Non-self-similar crack growth occurs be- 
cause of the anticlastic bending effect. The tangent stiffness 
matrix in the Newton-Raphson methods did not converge at 
the limit point because of the large value of the maximum 
interfacial strength, T c . The T c was reduced by half of its 
original value and a converged solution was obtained. Any of 
the modifications to the tangent stiffness matrix discussed in 
the section of interface elements, produced converged solu- 
tions without having to modify the originial value of T c . 

The ELS and ENF test specimens shown in Figure 5b and 
5c are used to determine the interlaminar' fracture toughness 
in Mode II. For the ELS, the load P is applied at the tip such 
that the lower arm of the ELS remains in contact with up- 
per arm. The other end of the specimen is clamped. The ENF 
specimen is simply supported, and the downward vertical dis- 
placement W 2 is specified at the mid-span of the specimen. 
The corresponding reaction force Pz is computed. The re- 
sponse of the ELS and ENF is shown in Figure 8a and 8b. For 
a loading-unloading cycle, excellent agreement of the FEM 
results are obtained compared to the closed form solutions. 

The FRMM test specimen is shown in Figure 5d, and is 
used to evaluate empirical failure criteria for mixed-mode de- 
lamination. The displacement w is specified at the tip of the 
upper arm and the corresponding reaction force P is com- 
puted. Mode I is 43% and Mode II is 57%. The response 
for a = 2 and a = 4 is shown in Figure 9a and 9b respec- 
tively. For a loading-unloading cycle, excellent agreement of 
the FEM results are obtained compared to the closed form 
solutions. 

The MMB test specimen is shown in Figure 5c, and is 
used to evaluate empirical failure criteria for mixed-mode de- 
lamination. The length of the lever arm c described in the 
report by Reeder 20 is chosen such that the mixed mode ra- 
tio from pure Mode I to pure Mode II can be varied. In 
this paper, c = 43.72 mm, so that the Mode I and Mode II 
contributions are 50% each. The MMB is simply supported, 
and two proportional loads are applied. The load Pi is ap- 
plied upward at the tip of the upper arm, and another load 
Pz is applied downward at the mid-span. During the load- 
ing, the ratio Pi/Pz = c/(c + L) is fixed. The responses for 
a = 4 and a = 2 are shown in Figure 10a and 10b. The fi- 
nite element response is compared to the analytical solutions 
in the appendix. In the first analysis, geometric nonlinear- 
ity is used. In the second analysis both geometric linearity 
and nonlinearity are compared with the analytical solutions. 
The discrepancies on the response corresponding to the stable 
crack growth of the load-deflection response are because the 
analytical solution does not consider the effects of geomet- 
ric nonlinearities. Excellent agreement is obtained with the 
analytical solutions. 


CONCLUSIONS 

An irreversible constitutive law that describes the delam- 
ination process is presented. The constitutive law is imple- 
mented with interface element to predict delamination. It pre- 
dicts initiation of delamination based on a multi-axial stress 
criteria, and progression of delamination based on an em- 
pirical fracture criteria. A damage parameter is included to 
prevent the restoration of the previous cohesive state between 
the interfacial surfaces. To demonstrate the irreversibility ca- 
pability of the constitutive law, steady-state crack growth is 
simulated for quasi-static loading-unloading cycle of various 
fracture test specimens. The finite element solutions are in 
excellent agreement with the analytical solutions. 
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APPENDIX 

The beam analytical solutions based on linear elastic frac- 
ture mechanics for the MMB test specimen are presented 
without details. In general, the total energy release rate is 

Gt = Gj + Gu (44) 

Gi and Gu are the Mode I and Mode II energy release rate 
contributions. The delamination propagates when, 

G t = G c = GT + GJj (45) 

and G c is the critical energy release rate, G™ and Gf] are 
the the Mode I and Mode II energy release rates at crack 
propagation. For all the fracture test specimens, it is possi- 
ble to express <j> = Gf /Gfj, where 4> € [0, oo), so that the 
G c value can be computed based on the fracture criterion in 
Equation (2) 


G c 


(! + </>) 



The derivations to obtain the expression of 4> for the MMB 
specimen are omitted here, and is 


Gj_ 

Gjj 


G 1 } 

7T 


a 


4 /6c- L 
3 V 2c+L 


(47) 


where c is the length of the lever arm, and L is the length of 
the MMB specimen. For simplifying purposes, the loads Pi 
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and Pu associated to modes I and II respectively are defined 
as 


p, = L - 

c 


6c — L 
4 L 


Pi 


Pn = - 

c 


2 c+ L 
L 


Pi (48) 


The load Pi is defined in Figure 5d. The initial load- 
deflection response is linear and given by 


W\ = 


16L 

37 


4 L 


i a o 

El 


(49) 


where E is the Young’s Modulus and I is the moment of in- 
ertia. The load-deflection response, when delamination prop- 
agates with a < L/2 is 


16P/ / 8 BEIG C \ 3/2 
3 El \64Pj + 3 Pfj J 


(50) 


where B is the width of the beam. The load-displecement 
relation when delamination propagates with a > L/2 is ob- 
tained by solving the quadratic equation for a. 


(64 Pf + 3 Pfj - 64 P/P//) a 2 (5 1 ) 

- (6Pf [ - 32 P 7 P 7 /L) a - (3 PfjL 2 - 8BEIG C ) = 0 


and substituting its solution into 


IV 1 


16L / 6 c — L\ Pi a 3 

77 V 4L J ~eT 


(52) 
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